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provided  by  NASA/Ames  Research  Center. 
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Scientific  Officer.  Dr.  Robert  S.  Rogallo  is  the  NASA  advisor. 
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INTRODUCTION 


Under  the  sponsorship  of  the  Office  of  Naval  Research  (ONR) , 
Nielsen  Engineering  & Research,  Inc.  (NEAR)  is  conducting  a program 
the  objective  of  which  is  the  testing  of  turbulence  models  using 
the  most  accurate  methods  of  computing  turbulent  flows  now  available. 
In  order  to  limit  the  number  of  possible  sources  of  discrepancy 
between  the  predictions  of  a model  and  the  results  of  an  accurate 
computation,  the  program  has  begun  by  looking  at  the  simplest  tur- 
bulent flows  and  the  simplest  models.  Since  the  program  is  expected 
to  continue  for  a few  years,  it  has  been  designed  so  that  new  fea- 
tures (geometric  complexity  and/or  physical  phenomena)  will  be 
added  one  at  a time.  This  should  provide  an  expanding  base  of 
confidence  to  build  on  and  will,  we  hope,  avoid  some  of  the  dif- 
ficulties that  other  researchers  have  had  in  sorting  out  various 
effects  in  turbulent  flows  and  in  learning  to  model  them. 

There  are  two  distinct  levels  of  turbulent  flow  computation  for 
which  models  will  be  investigated.  When  the  time-averaged  Navier- 
Stokes  equations  are  used  to  describe  a flow,  the  Reynolds  stresses, 
which  are  essentially  the  averages  of  products  of  the  fluctuating 
components  of  the  velocity,  need  to  be  modeled.  In  large-eddy 
simulation,  on  the  other  hand,  averages  of  the  products  of  the 
small-scale  components  occur  and  require  modeling.  By  analogy,  these 
are  called  the  subgrid-scale  Reynolds  stresses  but  they  represent  a 
different  set  of  physical  phenomena  than  the  time-average  Reynolds 
stresses.  However,  it  is  believed  that  both  sets  of  quantities  can 
be  modeled  in  similar  ways. 

The  basic  approach  to  model  validation  used  is  to  compute  an 
accurate  estimate  of  the  quantity  to  be  modeled  and,  simultaneously , 
the  value  that  the  model  would  predict  for  the  same  quantity.  Com- 
parison, usually  by  means  of  a correlation  coefficient,  then  pro- 
vides the  information  as  to  the  validity  of  the  model;  constants  in 
a model  can  also  be  evaluated.  This  general  approach  could  be 
implemented  in  several  ways.  In  the  first,  exact  Navier-Stokes 
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simulations  of  turbulent  flows  could  be  used  to  validate  both  subgrid- 
scale  and  time-average  models.  The  difficulty  is  that  only  a few 
turbulent  flows  (all  at  low  Reynolds  number)  are  accessible  to  exact 
simulation  and  this  severely  restricts  what  can  be  done  in  terms  of 
time-average  model  testing.  In  another  possible  implementation, 
large-eddy  simulations  could  be  used  to  test  time-average  models; 
the  problem  here  is  that  the  effects  of  the  subgrid-scale  turbulence 
(which  has  to  be  modeled  in  the  large-eddy  simulation)  may  be  dif- 
ficult to  estimate. 

Since  it  is  the  implementation  of  the  basic  approach  that  pro- 
vides the  best  accuracy,  we  have  chosen  to  begin  by  using  exact 
Navier-Stokes  simulations  to  test  subgrid-scale  models.  We  believe 
that  this  is  the  area  in  which  the  most  information  can  be  generated 
in  the  shortest  time.  One  can  also  argue  that  the  same  models  ought 
to  be  good  for  both  kinds  of  modeling  and  we  are  therefore  indirectly 
generating  information  about  time-average  modeling.  Of  course,  the 
last  statement  needs  to  be  checked  carefully  and  this  will  be  done 
in  later  stages  of  the  program. 

The  exact  simulations  that  we  are  using  as  the  basis  for  the 
model  checking  are  provided  by  Dr.  Robert  Rogallo  of  NASA/Ames 
Research  Center.  NASA  has  also  made  their  computer  available  to 
this  program  at  no  cost.  To  date,  the  only  flow  field  that  has  been 
available  to  us  has  been  a simulation  of  the  decay  of  homogeneous 
isotropic  turbulence.  We  have  used  these  data  as  input  to  a computer 
code  that  we  have  written  to  process  the  data  and  do  the  necessary 
correlations.  As  an  initial  test  of  both  the  input  data  and  our 
code,  we  compare  our  results  with  those  of  Clark,  et  al.  (ref.  1). 

The  results  of  these  tests  are  given  in  this  report  and,  although  a 
few  further  checks  are  necessary,  they  have  shown  that  both  Rogallo' s 
code  and  ours  seem  to  be  working  satisfactorily.  Extensions  of  this 
work  are  currently  underway. 

The  next  section  of  this  report  contains  a description  of  the 
computer  code  developed  by  NEAR.  The  third  section  presents  the 


6 


w 


results  we  have  generated  to  date,  and  the  last  section  describes 
extensions  of  this  work  to  be  undertaken  in  succeeding  contract 
years . 
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THE  NEAR  COMPUTER  CODE 


Description 


The  spectral  simulation  of  the  Navier-Stokes  equations  developed 
by  R.  S.  Rogallo  (ref.  2)  on  the  ILLIAC  IV  at  the  NASA/Ames  Research 

Center  results  in  values  for  the  three  velocity  components  u.  at 

3 1 

each  point  in  a 64  grid  at  each  time  step  in  the  calculated  evolu- 
tion of  a homogeneous  incompressible  flow.  If  the  Tay lor-microscale 
Reynolds  number  (R, ) of  the  flow  being  simulated  is  less  than  about 
40,  this  grid  is  fine  enough  to  capture  essentially  all  of  the  energy 
and  dissipation  in  the  flow,  and  these  results  can  be  considered 
"exact" . 


Imagine  that  the  same  homogeneous  incompressible  flow  is  to  be 
alculated  on  a coarser  grid  using  a large-eddy  simulation.  The 
equation  set  to  be  solved  is 


9u . 
l 

9x^ 


9u . 

l 

9t 


= 0 


(1) 


9 i 9P  , „2- 

-r (u.u.)  = - - + vV  U . - 

9x . l 3 9x . l 


1 " i J 

where  the  overbar  denotes  a spatially  filtered  variable 


(x • • + C.  . ) , 

3x,  ij  lj'j 


f (x) 


G(x-x')f (x')dx', 


(2) 


G (x)  is  the  selected  spatial  filter  function,  f (x)  is  any  field 
variable,  and  the  indicated  integration  is  over  all  space.  The 
underscore  denotes  a vector  quantity.  Also  in  equation  set  (1), 


P = p/p  + 1/3  (u£u£  + 2uRu^) 

Vj = “I3!  ■ 1/3  “K 

and  q j - GiU'  + ulGj  - 2/3  l^u') 
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The  solution  variables  in  these  equations  are  u^  and  P;  the  terms 

involving  finer  scales  than  are  resolvable  on  the  coarse  mesh  of 

the  LES  (i.e.,  those  involving  x^  and  must  be  modeled  in 

terms  of  u . . 

1 

The  object  of  our  work,  as  previously  stated,  is  to  evaluate 
models  of  the  type  used  in  large-eddy  simulations.  Our  evaluation 
is  of  several  eddy-viscosity  models  for  t. . (for  a model  for  C. . 

i]  ij 

and  its  evaluation,  see  ref.  1) . Because  in  large-eddy  simulations 
models  constructed  for  are  often  used  to  approximate  the  com- 

bination x. . + C. .,  we  also  evaluate  how  accurately  these  same 
il  il 

models  represent  this  quantity.  The  methodology  used  to  conduct 

this  evaluation  is  in  general  terms  as  follows.  Rogallo's  code  is 

3 

used  to  generate  the  velocity  field  u^  on  a 64  grid  (with  grid 
spacing  A)  at  a specified  time  step  in  the  evolution  of  a selected 
homogeneous,  incompressible  flow.  For  a selected  filter  function 
G(x)  (a  three-dimensional  Gaussian  with  filtering  length  scale 
A_  = 8A  is  the  initial  choice  in  this  work) , exact  values  are  calcu- 

a - 3 

lated  for  u.,  u!  and  x. . on  the  64  grid  using  relations  (2)  and  (3) 


and  this  u.  field. 
i 


il 

Values  of  lx. ./9x.  are  calculated  on  the  64' 

-3  il  1 


grid.  A coarse  grid  (16  , with  spacing  A = 4A,  thus  A /A  = 2) 

C cl  C 

representing  that  used  in  a hypothetical  LES  of  the  same  flow  is 

3 - 

overlaid  on  the  64  grid.  Exact  values  of  u.,  x. . and  9x. ./9x.  on 

3 i il  il  1 

the  coarse  grid  are  extracted  from  the  64 J grid.  Using  these  fields, 

exact  values  of  u.  9x. ./  x.  are  calculated  on  the  coarse  grid. 

i il  1 

Models  (M. .)  for  x. . in  terms  of  u.  are  also  calculated  on  the 
i]  il  i 

coarse  grid,  as  are  values  of  9M.  ./3x.  and  u.  9M . ,/9x..  The  deriv- 
y il  1 i il  1 

atives  required  for  the  model  calculations  are  done  using  any  one  of 

several  numerical  schemes.  At  this  point,  we  possess  (on  the  coarse 

grid)  both  exact  and  modeled  values  of  the  subgrid-scale  stress 

(x..  and  M..),  the  divergence  of  this  stress  (9x../9x.  and  3M../9x.), 
il  il  il  1 _ il  1 

and  the  energy  dissipation  by  this  stress  (u^  9x^/9Xj  and  3M^/9x.) 


For  the  case  where  x^j  + C^j  is  being  modeled,  substitute  x^j  + Cij 

for  x..  in  this  description, 
il 
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The  assessment  of  the  validity  of  the  models  used  can  now 

proceed  on  three  levels  as  was  pointed  out  by  Clark  (ref.  1) : 

(1)  the  tensor  level,  where  models  (Mr  ^)  are  compared  directly  to 

the  exact  subgrid-scale  stress  (t^);  (2)  the  vector  level,  where 

9M.  ./9x.  is  compared  to  3x.  ./5x.;  and  (3)  the  scalar  level,  where 
13  3 13  3 

u.  3M.  ./9x.  is  compared  to  u.  9t.  ./3x..  At  each  of  these  levels, 
i 13  3 i 13  3 

the  comparison  is  done  by  means  of  calculating  correlation  coeffi- 
cients 


C ( M , X ) 


<MX- 


? 1/2  „ 1/2 

<M  > <X  > 


(4) 


* 1 r A A 

where  < > = — - l ( ),  M = model  value,  X = exact  value.  The 

16  16 3 

magnitude  |C(M,X) ' will  vary  between  0 (if  M and  X are  totally 
unrelated)  to  1 (if  the  model  is  exact  to  w’ithin  a multiplicative 
constant).  Notice  that  at  the  tensor  level,  a correlation  coeffi- 
cient is  calculated  for  each  stress  component  (6  of  which  are 
independent) ; the  arithmetic  average  of  these  coefficients  gives 
a sense  of  the  correlation  for  the  whole  tensor.  Similarly,  at 
the  vector  level,  there  are  three  correlation  coefficients  (and 
their  average),  and  at  the  scalar  level,  one  correlation  coefficient. 

These  correlation  coefficients  are  independent  of  the  value 
used  for  the  constant  in  a given  model.  Values  for  the  constants 

★ 

can  be  determined,  however,  by  forcing  agreement  of  the  root  mean 
square  (RMS)  modeled  and  exact  values.  At  the  tensor  level,  a 
separate  value  of  the  constant  can  be  obtained  in  this  way  for  the 
diagonal  and  off-diagonal  components;  another  value  of  the  tensor- 
level  constant  can  be  derived  by  applying  this  RMS  criterion  to 
the  sum  of  squares  of  the  stress  components  at  each  grid  point. 
Similarly,  three  values  for  the  constant  (and  an  "overall"  value) 
are  available  at  the  vector  level,  and  one  at  the  scalar  level. 


On  the  16 


3 


mesh 
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The  models  evaluated  are  all  of  the  eddy-viscosity  type: 


M . . = a . . - 1/3  a,  ,6  . . 
13  19  kk 


a.  . = 2KS . . 
13  13 


S . . = 1/2 
13 


9u . 
l 

9x  . 


9u  . 


9x . 

l 


(5) 


J 


The  models  are 
1. 


2. 


3. 


2 1/2 

The  Smagorinsky  Model,  K = (C  A ) [2S. .S. .]  ' , 

3 1 s a.  13  19  ‘ 


where  A. 


is  the  averaging  length  scale,  and  Cs  is  the  model  con- 
stant . 


The  Vorticity  Model,  K = (C  A ) ' (w  . oj  . ) 1//2 

v a ii 

the  model  constant  and  w . = e.  9u, /9x.. 

l 19k  k'  9 


where  C is 


1/2 


The  Kinetic  Energy  Model,  K = (C  A /3)  (u'u,') 

-1  g a'  k k 

C is  the  model  constant  and  1/3  u’u,'  is  the  exact 
q k k 

subgrid-scale  kinetic  energy. 


where 


4. 


The  constant-eddy-viscosity  model,  K = Cc . 


A program  has  been  written  for  the  CDC  7600  computer  to 
accomplish  the  above  tasks.  Because  even  this  computer  cannot 
store  in  core  a complete  velocity  field  on  a 64  grid,  this  program 
incorporates  a considerable  amount  of  data  transfer  to  and  from 
auxiliary  disc  memory.  This  data  transfer  results  in  lengthy  res- 
idence time  in  computer  input/output  queues,  although  computation 
times  are  reasonable.  Unfortunately,  it  is  usual  for  any  sophis- 
ticated computer  system  to  "craah"  quite  frequently,  and  it  so 
happens  that  the  mean  time  between  crashes  for  the  CDC  7600  is  less 
than  the  residence  time  to  complete  the  analysis  of  a time  step  for 
a given  flow  field.  It  was  necessary  therefore  to  divide  the  pro- 
gram up  into  ten  subprograms,  each  with  restart  capability  in  the 
event  of  a system  failure  during  execution.  These  subprograms  are 
run  sequentially  and  are  named  CALCIA,  CALCIB,  CALCIC,  CALCII, 
CALCIII , CALCIV,  CALCV , CALCVI , CALCVII  and  INCORE.  A description 
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of  these  subprograms  and  the  major  subroutines  used  can  be  found 
in  the  Appendix.  This  section  is  concluded  with  a brief  discussion 
of  the  means  used  to  verify  the  accuracy  of  the  results  of  these 
programs,  and  a discussion  of  the  changes  made  to  evaluate  the 
same  models  as  applied  to  + C^j. 


Checkout  of  the  Computer  Code 


To  aeuug  these  programs  and  verify  their  accuracy,  they  wer-e 
used  to  analyze  an  analytically  specified,  artificial  "flow  field", 
the  Taylor-Green  flow  field: 


u^  = cos  ax^  sin  ax2  cos  ax^ 
U2  = -sin  ax^  cos  ax2  cos  ax^ 


> 

J 


(6) 


where  a = 2tt/64,  x^  are  in  the  domain  [0,63].  Notice  that  this 

"flow  field"  satisfies  continuity  exactly,  but  does  not,  of  course, 

have  any  features  of  turbulence.  Its  utility  for  checking  out  the 

code  lies  in  its  simplicity,  its  exact  periodicity,  and  the  fact 

that  because  of  its  wave  number  content,  the  quantities  u^,  , 

3t. ./  x.,  S.,  F.,  S.  and  F.  can  be  calculated  exactly  by  the 
ij  j 111  1 

CDC  7600  program  and  compared  to  simple  closed-form  expressions. 

Quantities  used  in  the  model  calculations,  although  not  calculated 

exactly  on  the  16^  mesh,  may  also  be  compared  to  (more  complicated) 

closed-form  expressions  to  ensure  that  the  code  is  operating 

correctly.  These  comparisons  have  been  carried  out  on  randomly 

3 3 

selected  elements  of  the  appropriate  mesh  (64  or  16  ) . 

Notice  that  because  = 0,  = T3^  = T23  = T32  = U3  ^T3j^xj 

Correct  calculation  of  these  zero  elements  does  not  guarantee  that 
some  problem  is  not  being  masked.  Therefore,  the  Taylor-Green  flow 
was  analyzed  twice  more,  each  time  after  rotating  the  coordinate 
system  to  place  the  zeros  in  different  elements. 


0. 
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Modifications  to  the  Code  to  Apply  the  Same  Models 


To  i . . + 
ij 


This  was  simply  accomplished  by  changing  those  portions  of 
the  code  that  do  input/output  on  the  64^  grid  to  include  the 
appropriate  elements  of  C^j.  A second  version  of  the  program 
set  was  prepared  in  which  these  changes  were  made. 


F 

I 

5: 


ki 
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RESULTS  TO  DATE 


Since  the  approach  to  model  validation  adopted  by  this  program 
requires  the  use  of  two  large  computer  codes  (one  to  calculate  the 
field  and  a second  to  analyze  the  data  and  produce  correlations) , it 
was  decided  early  on  that  the  initial  task  should  be  the  repetition, 
as  closely  as  possible,  of  at  least  one  case  that  has  already  been 
computed  and  analyzed.  The  only  case  available  at  the  present  is 
that  of  Clark,  et  al.  (ref.  1).  Clark  simulated  the  2.54  cm  grid- 
turbulence  experiment  of  Comte-Bellot  and  Corrsin  (ref.  3)  and,  on 
obtaining  good  agreement  with  the  experimental  data,  he  used  the 
computed  results  to  test  subgrid-scale  models.  We  will  present 
his  results  alongside  ours  below. 

To  facilitate  the  comparison,  we  wrote  the  computer  code  de- 
scribed in  the  previous  section  and  requested  Dr.  Rogallo  to  run  a 
simulation  of  the  Comte-Bellot  and  Corrsin  flow.  His  results  were 
provided  to  us  in  the  form  of  a tape  which  was  used  as  the  input  to 
our  program.  Before  going  on  to  the  results,  it  is  important  to 

compare  Rogallo' s simulation  with  Clark's.  Both  programs  compute 

3 

the  flow  field  on  a 64  grid.  Clark's  program  used  a fourth-order 
spatial  finite-difference  method  while  Rogallo 's  uses  Fourier 
methods  (in  fact  he  treats  the  Fourier  transform  of  the  velocity 
as  the  basic  dependent  variable  but  that  difference  is  not  impor- 
tant) . Clark  used  a third-order  predictor-corrector  scheme  for  the 
time  advancement  while  Rogallo  has  chosen  the  fourth-order  Runge- 
Kutta  method.  These  differences  are  not  expected  to  have  a sig- 
nificant effect  on  the  results.  For  ease  of  comparison  in  this 
report,  time  is  expressed  in  terms  of  Clark's  time  steps;  i.e., 
number  of  time  steps  = (time  in  seconds) /. 007 3 . 

The  initial  field  for  Rogallo 's  program  was  constructed  in  the 
same  way  as  Clark's.  However,  in  constructing  the  spectrum,  Rogallo 
has  used  different  subdivisions  in  wave  number  space  than  did  Clark. 
As  a result,  Rogallo' s spectra  look  lumpy  even  though  they  are  in 
fact  the  same  as  Clark's,  cf..  Fig.  1.  A more  serious  difficulty 
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is  that,  due  to  an  oversight,  a viscosity  was  used  that  is  twice  as 
big  as  it  ought  to  be  and  the  Reynolds  number  is  only  about  half  of 
the  desired  value  (see  Fig.  2).  Also,  the  time  scale  was  incorrect 
for  a similar  reason.  Consequently,  the  results  obtained  do  not 
match  the  experimental  energy  decay  or  dissipation  rate  as  well  as 
they  should;  these  results  are  shown  in  Figs.  3 and  4.  Note  that 
because  of  the  incorrect  time  scaling,  Rogallo's  last  time  point  is 
not  the  same  as  Clark's. 

We  now  turn  to  results  that  were  generated  by  our  code.  First 
we  give  some  of  the  statistics  of  the  flow  field.  The  skewness  of 
the  velocity  derivative  is  shown  in  Fig.  5 and  the  flatness  is  given 
in  Fig.  6.  In  Fig.  5,  we  compare  the  skewness  values  obtained  from 
each  of  the  three  velocity  components  with  those  of  Clark.  Clark 
provided  the  skewness  for  only  one  component  of  the  velocity  but 
has  informed  us  privately  that  the  other  components  show  about  the 
same  scatter  as  the  results  we  obtained.  Our  skewnesses  tend  to  be 
a bit  higher  than  his.  This  is  probably  due  to  the  use  of  Fourier 
methods  to  obtain  the  velocity  derivatives.  Fourier  methods  differ- 
entiate more  accurately,  especially  at  high  wave  numbers,  and  the 
skewness  is  sensitive  to  the  high -wave-number  components  of  the 
velocity  field.  Similar  remarks  apply  to  the  flatness  but  Clark  did 
not  report  values  of  this  quantity.  We  also  note  that  skewness  seems 
to  rise  to  the  equilibrium  value  more  slowly  in  our  calculation  than 
in  Clark's.  The  precise  reason  for  this  is  not  known  but  we  would 
guess  that  it  may  be  due  to  the  high-wave-number  portion  of  the  field 
requiring  more  time  to  come  to  equilibrium  than  the  low-wave-number 
part.  Since  our  calculation  of  the  skewness  is  more  sensitive  to 
the  high-wave-number  components,  this  might  explain  the  relatively 
slow  rise  of  the  skewness.  Also  shown  in  Figs.  5 and  6 are  the 
skewness  and  flatness  of  the  filtered  field.  As  expected,  they  are 
quite  a bit  below  the  values  for  the  full  field.  We  intend  to 
further  investigate  the  effect  of  filtering  on  the  skewness  and 
flatness  in  the  near  future.  This  would  assist  in  the  comparison 
of  large-eddy  simulations  with  experiment,  cf.,  Ferziger,  et  al. 

(ref.  4) . A preliminary  look  into  this  has  not  provided  anything 
useful  in  this  regard. 
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We  now  turn  to  results  obtained  from  the  filtered  field.  As 


explained  in  the  previous  section,  Rogallo's  velocity  field  was 
filtered  by  the  use  of  Fourier  methods  using  a Gaussian  filter. 

Once  the  filtered  field  has  been  computed,  we  can  obtain  the  subgrid 
scale  field  by  subtraction.  From  this  point,  we  can  go  on  to  cal- 
culate the  subgrid-scale  Reynolds  stresses,  the  derivatives  of  the 
filtered  field  and  all  of  the  quantities  that  go  into  a model.  Vie 
have  so  far  investigated  only  those  models  that  were  considered  by 
Clark:  the  Smagorinsky  model,  the  vorticity  model,  a one-equation 

subgrid-scale  kinetic  energy  model,  and  a constant-eddy-viscosity 
model . 

One  measure  of  the  overall  accuracy  of  the  calculations  comes 
from  comparing  dissipation  caused  by  the  subgrid-scale  Reynolds 
stresses.  These  are  displayed  in  Fig.  7.  The  comparison  with 
Clark's  results  is  fairly  good.  The  reason  for  the  higher  values 
at  the  early  time  steps  in  our  calculation  is  not  known  and  seems 
a little  surprising  in  view  of  the  fact  that  our  field  seems  to 
come  to  equilibrium  more  slowly  than  Clark's. 

Now  we  come  to  the  detailed  correlation  results.  As  mentioned 
earlier,  correlations  can  be  done  at  three  different  levels.  At  the 
most  detailed  level,  the  subgrid-scale  Reynolds  stresses  are  com- 
pared with  the  model  directly;  Clark  called  this  the  tensor  level. 

At  the  second  level  we  note  that  it  is  only  the  divergence  of  the 
stress  tensor  that  appears  in  the  dynamical  equations  and  we  compare 
the  divergence  of  the  exact  Reynolds  stress  with  the  divergence  of 
the  model;  Clark  termed  this  the  vector  level  of  comparison.  At 
the  crudest  level  we  argue  that  the  main  function  of  the  subgrid- 
scale  model  is  to  dissipate  the  kinetic  energy  of  turbulence  in  the 
right  amount  and  at  the  right  places  and  we  compare  the  energy 
dissipation  of  the  exact  result  and  the  model  prediction.  This  is 
done  by  taking  the  scalar  product  of  the  divergence  with  the  vel- 
ocity vector  itself;  Clark  calls  this  the  scalar  level  of  comparison 
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We  thus  have  four  different  models  that  we  wish  to  evaluate  at 
three  different  levels.  The  comparison  should  also  be  made  at 
several  different  time  steps  to  check  that  the  results  are  not 
peculiar  to  one  particular  field.  (They  ought  to  be  checked  on 
different  realizations  of  the  flow  as  well  but  the  expense  of  gen- 
erating flow  fields  has  precluded  this.  A small  amount  of  this 
kind  of  checking  will  be  done  in  the  future.)  The  question  of 
whether  there  is  a variation  with  time  is  easily  disposed  of. 

Figure  8 shows  the  results  for  the  Smagorinsky  model  using  a standard 
centered  difference  scheme  to  do  the  differentiation  required  in  the 
model  calculations,  and  it  is  seen  that  the  variation  with  time  is 
not  significant.  The  results  for  the  vorticity  and  constant-eddy- 
viscosity  models  are  very  similar  and  are  not  shown  for  this  reason; 
th‘£  latter  is  a bit  of  a surprise.  For  the  turbulent  kinetic  energy 
model,  however,  we  find  that  the  correlation  coefficient  is  essen- 
tially constant  but  there  is  a small  increase  in  the  constant  with 
time  as  shown  in  Fig.  9.  This  may  have  important  implications  for 
turbulence  modeling  but  the  result  was  only  recently  obtained  and 
we  have  not  had  time  to  analyze  it  as  yet.  In  view  of  the  constancy 
of  most  of  the  results  with  time,  we  will  concentrate  on  one  time 
step  from  now  on. 

In  Tables  1-4  the  correlation  coefficients  are  given  for  the 
four  models  at  each  of  the  three  levels.  There  is  no  significant 
difference  among  the  models  as  far  as  the  correlation  coefficients 
are  concerned.  In  fact,  at  the  tensor  level  the  correlations  are 
much  closer  when  a given  component  is  compared  for  the  various 
models  than  when  the  individual  components  are  compared  for  a par- 
ticular model.  This  is  consistent  with  Clark's  results  and  probably 
means  that  the  correlation  is  a property  of  eddy-viscosity  models 
in  general  and  has  little  to  do  with  the  particular  form  chosen  for 
the  eddy  viscosity.  If  this  is  so  (and,  again,  it  is  in  agreement 
with  Clark's  results),  it  means  that  the  prognosis  for  one-  or  two- 
equation  subgrid-scale  models  is  not  good--a  constant-eddy-viscosity 
model  will  do  essentially  as  well  as  a more  sophisticated  model. 
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Or  to  put  it  another  way,  a more  sophisticated  model  will  do  as 
badly  as  the  constant-eddy-viscosity  model.  This  result,  important 
though  it  may  be,  requires  further  checking  on  other  flows  before 
it  can  be  accepted  as  established.  The  generalization  of  this 
result  Reynolds-stress  modeling  is  on  even  shakier  ground  but 
it  sug je  s an  important  line  of  research  that  we  intend  to  take 
as  soon  is  the  tools  and  the  necessary  data  are  available.  The 
constants  for  each  of  the  first  three  models  (obtained  by  equating 
the  RMS  value  of  the  model  with  the  exact  result)  are  given  in 
Table  5. 

In  Table  6 are  presented  the  correlation  coefficients  obtained 
using  different  methods  of  calculating  the  derivatives  in  the  model 
calculations.  The  methods  used  are  the  standard  centered  difference 
formula  (L  = 1),  a centered  difference  formula  using  more  widely 
spaced  points  (L  = 2) , and  the  spectral  method.  It  is  seen  that 
the  correlation  coefficients  obtained  from  the  standard  centered 
difference  formula  and  from  the  spectral  method  are  essentially  the 
same,  while  those  from  the  alternate  centered  difference  formula 
are  lower.  Use  of  a centered  difference  formula  is  approximately 
equivalent  to  using  a "box"  filter  of  width  equal  to  the  distance 
between  the  points  at  which  the  function  is  evaluated.  Thus  the 
use  of  difference  formulas  with  different  mesh  spacing  is  a quick 
and  easy  (although  non-exact)  way  to  search  for  effects  of  averaging 
in  the  model  calculations.  It  has  been  suggested  by  Leslie  and 
coworkers  (ref.  5)  that  subgrid-scale  models  are  improved  by  aver- 
aging them  over  small  regions  of  space;  they  suggested  averaging 
over  a region  of  characteristic  length  of  about  twice  the  width  of 
the  filter.  Recent  work  of  the  Stanford  group  provides  indirect 
confirmation  of  this.  However,  the  results  in  Table  6 for  L = 2 
vs.  L = 1 indicate  that  the  correlation  is  not  improved  by  averaging 
the  model  over  more  than  one  filter  width.  A more  careful  study  of 
this  will  be  undertaken  in  the  near  future. 
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Table  6 also  compares  our  results  with  those  of  Clark.  Clark 
used  a fourth-order  finite  difference  formula  on  his  8J  coarse 
grid  to  do  the  model  calculations.  It  is  clear  that  our  correlations 
(excluding  those  using  L = 2)  agree  well  with  his  at  the  tensor  level 
but  are  smaller  at  the  vector  and  scalar  levels  of  comparison.  In 
fact,  the  improvement  that  Clark  found  in  going  from  the  tensor  to 
the  vector  level  seems  to  be  absent  in  our  results.  The  reason  for 
this  is  not  understood  at  the  present  time  but  some  of  the  studies 
to  be  made  in  the  near  future  should  shed  light  on  this.  That  there 
are  several  possible  causes  of  the  discrepancies  is  clear  upon 
examination  of  Table  7,  which  summarizes  the  parameters  and  methods 
used  in  the  present  work  with  those  of  Clark.  Those  differences 
thought  to  be  important  are: 

1.  Our  calculations  were  made  at  a Reynolds  number  that  is  only 
about  half  of  Clark's.  It  is  possible  that  we  are  seeing  a Reynolds 
number  effect.  Since  the  effect  of  Reynolds  number  is  one  of  the 
highest  priority  items  for  the  near  future,  this  possibility  will 

be  checked  out  soon. 

2.  Rogallo's  calculation  used  Fourier  methods  whereas  Clark's 
used  finite  differences.  It  is  possible  that  the  two  flow  fields 
differ  in  the  high  wave  number  part  of  the  velocity  field.  Since 
these  contain  most  of  the  contribution  to  the  subgrid-scale  turbu- 
lence component,  it  is  possible  that  the  source  of  the  difference 
lies  here. 

3.  We  used  a Gaussian  filter  and  Clark  used  a box  filter.  The 
effect  of  filter  type  on  the  results  is  a relatively  easy  item  to 
check  using  the  code  we  have  developed.  We  will  be  looking  not  only 
at  the  two  filters  mentioned  above  (the  Gaussian  and  the  box)  but 
also  at  a filter  that  corresponds  to  a sharp  cutoff  in  Fourier  space. 

4.  We  used  a finer  "coarse"  mesh  than  did  Clark  (16^  vs.  8 3) 

to  improve  the  accuracy  of  the  numerical  techniques  used  in  the  model 
calculations  and  to  increase  the  sample  size  for  our  statistical 
evaluation  process.  After  the  other  effects  listed  above  are  eval- 
uated, this  difference  can  be  eliminated  if  it  is  warranted. 
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Table  8 compares  our  model  constants  obtained  using  the  standard 
centered  difference  formula  and  the  spectral  method  with  those  that 
Clark  obtained.  The  model  constants  we  calculated  show  more  depen- 
dence on  the  method  of  differentiation  used  than  the  correlation 
coefficients  did.  With  the  exception  of  the  kinetic-energy  model, 
our  values  are  lower  than  Clark's.  We  suspect  that  this  may  be 
due  to  the  effects  listed  above,  but  it  is  impossible  to  be  con- 
fident about  this  until  the  further  checks  mentioned  above  are  made. 

Finally,  we  note  that  although  Clark  did  the  kind  of  testing 

done  above,  some  authors  define  the  subgrid-scale  Reynolds  stress 

to  be  u!u!  + u.u!  + u!u.  rather  than  u!u'.  We  have  therefore  tested 
il  il  il  il 

the  models  given  above  as  models  of  this  modified  subgrid-scale 

Reynolds  stress  and  the  results  are  given  in  Table  9.  It  is  clear 

that  the  models  work  better  for  the  u!u'.  term  than  for  the  combined 

i 1 

term.  Again,  this  result  may  depend  on  the  effects  listed  above. 
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FUTURE  DIRECTIONS 


As  stated  earlier,  this  work  is  only  the  beginning  of  a longer 
program.  The  computer  code  that  was  used  to  obtain  the  results  of 
the  previous  section  has  been  available  for  only  a short  time  and 
the  results  given  are  definitely  of  a preliminary  nature.  Essenti- 
ally, most  of  the  first  year  of  the  program  was  spent  in  developing 
a tool  which  will  be  used  in  the  future.  The  coming  year  should 
therefore  see  a considerable  increase  in  the  rate  at  which  the  work 
progresses.  A possible  bottleneck  may  be  the  fact  that  the  data  we 
analyze  are  obtained  on  the  ILLIAC  IV  by  other  researchers  and  is 
not  within  the  control  of  Nielsen  Engineering  and  Research,  Inc. 

Some  of  the  future  directions  were  already  stated.  We  will, 
in  the  near  future,  look  at  the  effects  mentioned  in  the  previous 
section.  Specifically,  the  effects  of  Reynolds  number,  filter  type 
and  filter  length  scale  will  be  looked  at.  If  necessary,  we  will 
also  look  at  a further  comparison  of  Clark's  and  Rogallo's  fields. 

The  work  will  also  be  extended  to  consider  flows  other  than 
decaying  homogeneous  turbulence.  Specif ically , we  will  look  at 
homogeneous  strained  and  sheared  turbulence.  We  will  also  begin 
to  look  into  the  possibilities  of  full  subgrid-scale  Reynolds  stress 
modeling,  i.e.,  the  possibility  of  treating  the  subgrid-scale  stiess 
by  means  of  a set  of  six  partial  differential  equations. 
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Time  Step, 

N[Time  = .0073  N sec] 


Figure  2.  Reynolds  number  based  on  Taylor  microscale  as 
a function  of  time,  Rogallo's  simulation. 


Dissipation  Rate 


Skewness 


Time  Step,  NlTime  = .0073  N seel 


Figure  5.  Skewness  as  a function  of  time. 
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Figure  7.  Subgrid-scale  dissipation  rate,  t . 
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Time  Step,  N[Time  .0073  N sec] 


Figure  8.  Time  development  of  correlation  coefficients  and  mode] 
constants,  Smagorinskv  model,  standard  centered  differences 
used  in  model  calculations. 


TABLE  1.  DETAILS  OF  CORRELATIONS,  SMAGORINSKY  MODEL,  TIME 
STEP  22.4,  STANDARD  CENTERED  DIFFERENCES  USED  IN 
MODEL  CALCULATIONS 


Tensor  Level 


\j-  1 

2 

3 

i = 1 

-.358 

-.345 

-.206 

2 

-.345 

-.288 

-.245 

3 

-.206 

-.245 

-.285 

Average  = -.280 


Vector  Level 

i = 1 2 3 

-.259  -.276  -.239 

Average  = -.258 

Scalar  Level 


TABLE  2.  DETAILS  OF  CORRELATIONS,  VORTICITY  MODEL,  TIME 
STEP  22.4,  STANDARD  CENTERED  DIFFERENCES  USED  IN 
MODEL  CALCULATIONS 


Tensor  Level 


\ 

j = 1 

2 

3 

i = 1 

-.326 

-.360 

-.200 

2 

-.360 

-.258 

-.228 

3 

-.200 

-.228 

-.260 

Average  = -.269 


Vector  Level 


1 2 3 

-.267  -.272  -.247 


Average  -.262 


Scalar  Level 
-.520 
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TABLE  3.  DETAILS  OF  CORRELATIONS,  SUBGRID-SCALE  KINETIC  ENERGY 
MODEL,  TIME  STEP  22.4,  STANDARD  CENTERED  DIFFERENCES 
USED  IN  MODEL  CALCULATIONS 


Tensor  Level 


\j = . 1 J 

2 

3 

i = 1 

-.379 

-.380 

-.230 

2 

-.  380 

-.319 

-.262 

3 

-.230 

-.262 

-.306 

Average  = -.305 


Vector  Level 


1 2 3 

-.285  -.304  -.260 


Average  = -.283 


TABLE  4. 
MODEL, 


DETAILS  OF  CORRELATIONS,  CONSTANT-EDDY-VISCOSITY 
TIME  STEP  22.4,  STANDARD  CENTERED  DIFFERENCES 
USED  IN  MODEL  CALCULATIONS 


Tensor  Level 


j = 1 

2 

3 

i = 1 

-.358 

-.362 

-.226 

2 

-.362 

-.301 

-.256 

3 

-.226 

-.256 

-.290 

Average  = -.293 


Vector  Level 


1 2 3 

-.270  -.297  -.251 


Average  = -.273 


Scalar  Level 


TABLE  5.  DETAILS  OF  MODEL  CONSTANTS,  TIME  STEP  22.4, 
STANDARD  CENTERED  DIFFERENCES  USED  IN  MODEL  CALCULATIONS 


Model 

Tensor  Level 

Smaqor insky 

Vorticity 

Subgrid-scale 
kinetic  energy 

Constant  for  diagonal 
elements 

.160 

.177 

.229 

Constant  for  off- 
diagonal  elements 

.167 

.183 

.245 

Constant  for  all 
elements 

.164 

.180 

.238 

Vector  Level 

Constant  for  i = 1 

. 189 

.209 

. 320 

Constant  for  i = 2 

.193 

.212 

. 324 

Constant  for  i = 3 

.193 

.213 

. 327 

Overall  constant 

.191 

.211 

.324 

Scalar  Level 

Constant 

.142 

.155 

.175 
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TABLE  6.  DEPENDENCE  OF  AVERAGE  CORRELATION  COEFFICIENTS  ON 
METHOD  OF  CALCULATING  DERIVATIVES  IN  MODELS, 

TIME  STEP  22.4 


Average  Correlation  Coefficients 

Centered  Difference  Spectral  Results  of 

Formula  Method  Clark  (ref.  1) 

(see  sketch  below)  


L = 1 L = 2 


Tensor  Level 


Smagorinsky 

-.280 

-.161  -.278 

-.277 

Vorticity 

-.269 

-.187  -.265 

-.260 

Kinetic  Energy 

-.305 

-.214  -.299 

-.303 

Constant  K 

-.293 

-.202  -.288 

-.295 

Vector  Level 

Smagorinsky 

-.258 

-.194  -.241 

-.346 

Vorticity 

-.262 

-.214  -.245 

-.327 

Kinetic  Energy 

-.283 

-.230  -.269 

-.362 

Constant  K 

-.273 

-.227  -.260 

-.356 

Scalar  Level 

Smagorinsky 

-.516 

-.474  -.492 

-.580 

Vorticity 

-.520 

-.496  -.496 

-.582 

Kinetic  Energy 

-.549 

-.506  -.537 

-.606 

Constant  K 

-.533 

-.502  -.523 

-.605 

3f(I) 

f ( I+L)  - f(I-L) 

3x 

2 (L) (Ax) 

1-2  I- 

1 I 1+1  1+2 

Centered 

Difference  Formula 
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TABLE  7. 


SUMMARY  OF  PARAMETERS  AND  METHODS  OF  PRESENT 
WORK  AND  THOSE  OF  CLARK  (REF.  1) 


Navier-Stokes  Solution 
Grid  (Spacing  = A) 
Space  Differencing 

Time  Differencing 

Initial  Energy 
Spectrum 

(RA ^ initial 

Filtered  Fields 
Grid 
Filter 

Filtering  length 
scale  (A  ) 

a 


Clark  (ref.  1) 


Present  Work 


64 


3 


64 


3 


fourth-order  spectral 

finite  difference 

third-order  fourth-order 

predictor-corrector  Runge-Kutta 

from  ref.  3 same  as  Clark's 


38.1 


22.3 


8 


3 


16 


3 


Box 


Gaussian 


8 A 


8A 


Model  Derivatives 


fourth-order 
finite  difference 


variable 
(see  text) 


TABLE  8.  DEPENDENCE  OF  MODEL  CONSTANTS  ON  METHOD  OF 
CALCULATING  DERIVATIVES  IN  MODEL,  TIME  STEP  22.4 


Model  Constant  (Tensor  and 
obtained  using 

Vector  values 
all  elements) 

Standard 

Centered  Difference 
Formula 

Spectral 

Method 

Results  of 
Clark  (ref.  1) 

Tensor  Level 

Smagor insky 

.164 

.128 

.247 

Vorticity 

.180 

.141 

.275 

Kinetic  Energy 

.238 

.186 

.175 

Vector  Level 

Smagorinsky 

.191 

.119 

.264 

Vorticity 

. 211 

.130 

. 247 

Kinetic  Energy 

. 324 

.163 

. 155 

Scalar  Level 

Smagorinsky 

. 142 

. 094 

.171 

Vorticity 

.155 

. 102 

.191 

Kinetic  Energy 

.175 

.100 

. 095 
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TABLE  9.  COMPARISON  OF  RESULTS  OBTAINED  BY  APPLYING  SAME 
MODELS  TO  Tji  + Cij  AND  TO  Tij,  TIME  STEP  22.4,  STANDARD 
CENTERED  DIFFERENCES  USED  IN  MODEL  CALCULATIONS 


Modeled  Quantity: 


u ! u ! 
i 3 

u ! u + u . u 
il  13 

+ u '•  u . 
i 3 

Tensor  Level 

Average 

Correlation 

Coefficient 

Model 

Constant 

Average 

Correlation 

Coefficient 

Model 

Constant 

Smagorinsky 

-.280 

.164 

-.198 

. 301 

Vorticity 

-.269 

. 180 

- .190 

. 332 

Kinetic  Energy 

-.305 

.238 

-.210 

.806 

Constant  K 

-.293 

— 

-.201 

— 

Vector  Level 

Smagorinsky 

-.258 

. 191 

-.178 

. 333 

Vorticity 

-.262 

.211 

-.187 

. 368 

Kinetic  Energy 

-.283 

. 324 

-.194 

.980 

Constant  K 

-.273 

— 

-.187 

— 

Scalar  Level 

Smagorinsky 

-.516 

.142 

-.483 

.216 

Vorticity 

-.520 

. 155 

-.493 

.236 

Kinetic  Energy 

-.549 

.175 

-.514 

.404 

Constant  K 

-.533 

— 

-.499 

— 
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APPENDIX  A 


SUBPROGRAM  AND  MAJOR  SUBROUTINES  IN  THE  NEAR  COMPUTER  CODE 

CONVERT 

Before  embarking  on  the  calculations  described  earlier,  the 

tape  generated  on  the  ILLIAC  (64  bits/word)  must  be  converted 

(on  the  CDC  7600)  to  the  CDC  7600  word  structure  (60  bits/word) 

and  separate  permanent  files  set  up  on  a private  disc  pack  for 

the  three  velocity  components  u..  These  steps  are  achieved  in 

program  CONVERT.  The  three  files  for  u.  each  contain  64  J (262,144) 

* ^ 
words.  A "header"  file  of  global  quantities  calculated  by  the 

ILLIAC  is  printed  out.  The  velocity  components  u^  are  in  the 

dimensionless  form  computed  on  the  ILLIAC  (the  normalization  is 

described  in  ref.  2). 

MAIN 

This  is  the  driver  program  for  the  CALC  series  of  subprograms. 

In  this  program,  certain  initializations  are  performed,  a back-up 
copy  of  the  current  version  of  file  LCMDAT  (described  below)  is 
created,  three  scratch  files  are  created  on  a disc  each  with  the 
capacity  for  a complex  64~*  field  (524,288  words),  and  control  is 
passed  to  the  appropriate  CALC  subprogram.  When  the  selected  CALC 
has  returned  control  to  MAIN,  transfer  of  the  results  of  that  CALC 
to  a new  version  of  file  LCMDAT  is  accomplished.  In  the  event 
of  system  failure  during  any  CALC,  the  back-up  copy  of  LCMDAT  can 
be  used  to  restart  that  CALC. 

CALCIA 

Permanent  files  on  the  private  disc  pack  are  set  up  to  later 

3 

accommodate  u.,  u!  and  -u.'u.’/l.  Each  of  these  files  is  64  words 
ii  k k 

long.  File  LCMDAT  is  also  established  on  the  private  disc  for 

* 

The  file  lengths  quoted  are  actually  slightly  larger  in  practice 
by  the  small  amount  required  for  system  overhead  (i.e.,  index  arrays). 
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later  use.  This  file  holds  values  on  the  16  mesh  of  9t  • /9x.,  u., 

^3  It  i 

u'u'  and  the  6 independent  components  of  t . . and  is  thus  13X16-5  = 
k k il 

53,248  words  long. 

In  this  CALC,  u^,  u^  and  the  skewness  and  flatness  for  u^  and 
U1  ^Sl'  Sl'  Fl'  Fl^  are  calcu-*-atec^  • These  quantities  are  defined 
as  follows: 


S . = - < 
i 


F.  = - < 

l 


2 3/2 


S . = 
l 


F . = 
l 


< T— - >/< 

r9u . i 4 
"(^1  V< 


3/2 


T-  2 2 

9u . , 
-^1  > 
9xi  J 


(A.l) 


In  these  definitions,  the  summation  convention  is  suppressed  and  < > 

3 

denotes  averaging  over  the  64  field.  The  skewness  and  flatness 
values  are  printed  out,  u^  and  uj!  are  stored  in  the  64^  files,  and 
appropriate  elements  of  u^  are  stored  in  LCMDAT . 

The  convolution  to  calculate  the  filtered  field  as  expressed  in 
eqn.  (2)  in  the  body  of  this  report  is  done  using  a system-provided 
Fast  Fourier  Transform  (FFT)  and  the  discrete  convolution  theorem 
(ref.  6) . The  required  multiplication  in  wave  space  is  done  in  sub- 
routine FILTER.  Differentiation  is  also  done  using  spectral  methods 
in  subroutine  DIFFER. 

The  calculations  in  CALCIA  (and  the  other  CALC’s)  are  done  in 
dimensionless  variables.  Throughout  this  series  of  programs,  veloc- 
ities remain  as  normalized  by  Rogallo,  and  space  variables  are 

3 

normalized  using  A,  the  grid  spacing  of  the  64  grid,  as  the  charac- 
teristic length. 

Operational  requirements  of  the  CDC  7600  dictate  that  the  calcu- 

2 

lations  in  all  CALC ' s be  done  in  a single  x-y  or  x-z  plane  (64 
elements) . Further,  the  FFT  can  only  be  invoked  with  respect  to  the 
dimensions  present  in  the  plane  currently  in  core  (i.e.,  if  an  x-y 
plane  is  in  core,  only  the  x and  y transforms  of  that  plane  can  be 
taken).  Thus,  multiple  passes  through  the  data  files  are  required. 
Input-output  routines  GETPL , GETRPL , PUTPL,  and  PUTRPL  do  most  of 
the  required  data  transfer. 


CALC IB 


In  a similar  fashion  to  that  described  in  CALCIA,  u 2,  and 

S^,  F 


F„  are  calculated.  The  skewness  and  flatness  values 

3 

u0  and  u^  are  stored  in  64  files,  and  appropriate 


S2 ' 13  2 ' l'2'  "2 
are  printed  out 

elements  of  u 2 are  stored  in  LCMDAT . 
CALCIC 


As  above,  u^,  u^,  S3,  S3,  F^ , F^  are  calculated  and  stored  or 
printed  out. 


CALC I I 


Values  of  u^u^,  3u^U2/3x1,  and  3u^u2/3x2  are  calculated  on  the  b4 

grid.  Appropriate  elements  (on  the  16J  grid)  are  stored  as  and  a. 

partial  values  of  3x~./3x.  and  3x.  ./3x.,  respectively,  in  LCMDAT. 

^ 3 3 4 3 3 

CALC I I I 


Values  of  U3U3,  Su^u^/Sx^,  and  Su^u^/Sx^  are  calculated  on  the 


3Tl3/3xi  is  stored  as  a 


64  grid.  In  LCMDAT,  u|u3  is  stored  as  X33 

partial  value  of  3x3j/3x^,  and  8X33/3X3  is  added  to  the  previous 

partial  value  of  3x. ./3x.. 

13  3 

CALCIV 


Values  of  u^u^,  Su^u^/Sx.,,  and  3u2U3/3x3  are  calculated  on  the  64 
grid.  In  LCMDAT,  ul,u^  is  stored  as  x23,  3x23/3x2  is  added  to  9l3j^,Xj' 
and  8X23/8X3  is  added  to  3x2j/3x^.. 


CALCV 


On  the  64  grid,  -1/3  u^u^,  u^u^,  “ 1/3  UkUk  = T11‘ 


and 


SXif/SXi  are  calculated.  The  first  field  is  stored  in  the  file 
defined  for  it  in  CALCIA.  In  LCMDAT,  133  and  u^~u^  are  stored,  and 
8X33/3X3,  is  added  to  8X3 j / 3 x ^ . 


CALCVI 


3 2 

On  the  64  grid,  u2 


1/3  u/u,'  = x„_  and  3x_0/3x0  are  calcu- 
k k 22  22  2 


lated.  In  LCMDAT,  x22  is  stored  and  3x22/3x2  is  added  to  3x2-/3xj1 
CALCVI I 

On  the  64  3 grid,  u'^2  - 1/3  u£u^  = x33  and  8X33/8X3  are  calcu- 
lated. In  LCMDAT,  X33  is  stored  and  8X33/8X3  is  added  to  3x3j/3x  . 
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INCORE 


LCMDAT  is  the  input  to  this  subprogram  which  operates  only 
3 

on  data  on  the  16  grid.  The  correlation  coefficients  and  model 
constants  discussed  earlier  are  calculated.  The  derivatives 
required  in  the  model  calculations,  for  example  to  calculate  . 
or  3M^j/9Xj,  are  done  by  the  spectral  method,  by  a standard 
centered  difference  scheme  (L  = 1 in  the  sketch  below) , or  by  a 
modified  central  difference  scheme  with  L > 1. 


-►  X 


1-2  1-1 


1+1  1+2 


Dimensional : 


9_f  1 _ f ( I+L)  - f(I-L) 

3xJ  ' 2^>Ac 


Dimensionless : 


f ( I+L)  - f(I-L) 
2(L)Ac/A 


J 


Ac  is  the  grid  spacing  on  the  coarse  (163)  mesh  (A  /A  = 4) 


(A. 2) 


FILTER 

This  subroutine  does  the  multiplication  in  wave  space  cor- 
responding to  a convolution  in  physical  space.  It  forms  the 
product 

f ( i , j ) = f(i,j)  g-^i)  g2(j)  g3(k)/643  (A. 3) 

on  the  kth  x-y  wave  number  plane,  or  the  product 

f (i , k)  = f ( i , k)  gi(i)  g2(j)  g3(k)/643  (A. 4) 

on  the  jth  x-z  wave  number  plane,  where  i,j,k  are  the  indices  in 
the  x,y,z  wave  number  directions,  respectively,  (')  denotes  a 
three-dimensionally  Fourier  transformed  variable,  and  g^,  g2»  g3 

4 3 


are  the  one-dimensional  Fourier  transforms  of  the  components  of  the 
filter  function,  G^ (x) , where 

G (x)  = G1(x1)  G2(x2)  G3(x3)  . (A. 5) 

This  construction  allows  for  the  possibility  of  a nonisotropic 
filter.  g^ , g2,  and  g^  are  generated  in  subroutine  GHATGEN. 


GHATGEN 

This  subroutine  forms  the  one-dimensional  Fourier  transforms 
of  the  filter  function  components.  The  filter  function  implemented 
at  this  time  is  a three-dimensional  Gaussian: 


G (x)  = 


/6A 


Aa/Aj 


exp  [- 


(A  /A) 

d 


, 2 . 2 
(X1  + X2 


X3> 


(A. 6) 


This  filter  is  isotropic  with  a filtering  length  scale  of  A . The 
value  of  A /A  is  taken  to  be  8,  implying  that  A /A  = 2.  The  set 
g.  is  the  set  of  discretized  continuous  one-dimensional  Fourier 
transforms  of  this  filter  function: 


gx  (i) 
g2(i} 

g3  (k) 


exp  [- 


exp  [- 


exp  [- 


( A / A) 

a 

ui 

6 

[64 

( A / A ) 2 

a 

6 

’ui' 

[64 

(A  /A)2 

cl 

7ik' 

6 

[64.1 

f 


J 


(A. 7) 


In  general,  the  functions  g^  are  constrained  to  be  real,  which 
implies  that  the  filter  function  must  be  even. 


DIFFER 

This  subroutine  performs  the  multiplication  in  wave  space 
corresponding  to  differentiation  in  physical  space.  It  forms  the 
product 
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l«i-5>  ■ ^ 


2rti 


64 


/\ 

f(i,j)  or  |p-(i,  j) 


*/-l  ^4-  f (i  / j ) (A.  8) 


in  an  x-y  wave  number  plane, 
/\ 


fl'i.k)  • /=! 


2ui 


64 


/\ 

f (i, k)  or  ||-(i,k) 


/=T  f(i,k)  (A. 9) 


in  an  x-z  wave  number  plane.  Note  that  the  indicated  multiplica- 
tions are  complex. 


GETPL  and  PUTPL 

These  subroutines  are  input  and  output  subroutines,  respec- 
tively, which  move  planes  of  complex  elements  between  small  core 
memory  and  disc  storage  using  random  access  methods.  The  file 
and  array  names  involved  and  the  plane  number  are  arguments  in  the 
call.  The  plane  can  be  either  an  x-y  plane  or  an  x-z  plane. 


GETRPL  and  PUTRPL 

These  subroutines  are  input  and  output  subroutines,  respec- 
tively, which  move  plane  of  elements  between  small  core  memory  and 
disc  storage,  just  as  above.  The  difference  is  that  for  GETRPL 
and  PUTRPL  the  elements  in  mass  storage  are  real,  while  the  array 
in  core  may  be  complex.  On  input,  this  means  that  the  imaginary 
part  of  the  array  is  set  to  0;  on  output,  the  imaginary  part  is  lost. 
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LIST  OF  SYMBOLS 


constant  in  constant-eddy-viscosity  model 
constant  in  subgrid-scale  kinetic-energy  model 
constant  in  Smagorinsky  model 
constant  in  vorticity  model 
cross-scale  stress,  eqn.  (3) 


3 2 

three-dimensional  energy  spectrum,  cm  /sec 


velocity-gradient  flatness,  eqn.  (A.l) 

velocity-gradient  flatness  for  filtered  field,  eqn.  (A.l) 
any  field  variable 
spatial  filter  function 

component  of  spatial  filter  function,  eqn.  (A. 5) 


one-dimensional  Fourier  transform  of  G^(x^),  eqn.  (A. 7) 


indices  in  x,  y,  and  z wave  number  directions 
eddy  viscosity 


shell  wave  number,  cm 


-1 


spacing  variable  in  centered  difference  formula,  eqn.  (A. 2) 

1 


subgrid-scale  stress  model,  M.  . = a.  . - T a.  . 6 . . 
^ lj  lj  3 kk  13 

reduced  pressure,  eqn.  (3) 


3u  . 
1 


3x_ 


pressure 

Reynolds  number  based  on  Taylor  microscale 

rate  of  strain  tensor  for  filtered  field,  S . . = i 

1 j z 

velocity-gradient  skewness,  eqn. (A.l) 

velocity-gradient  skewness  for  filtered  field,  eqn.  (A.l) 
time 

velocity 


<u 


Dx  . 

1 
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LIST  OF  SYMBOLS  (concluded) 


filtered  velocity 

subgrid-scale  velocity,  eqn.  (3) 

spatial  coordinate 

a.  . = 2KS.  . 

ID  ID 

grid  spacing  for  64  grid 

grid  spacing  for  coarse  grid 

length  scale  for  filter  function  .. 

x NJ_  9t  . . 

subgrid-scale  dissipation,  ex  = — j £ u.  ^ 1 ^ 

N 1 xj 

N = number  of  grid  points  in  coarse  grid 
kinematic  viscosity 
density 

subgrid-scale  stress,  eqn.  (3) 
vorticity 


Subscripts 


1,2,3^ 
i / j / k y 


coordinate  directions 


vector  quantity 


Superscripts 


L 


filtered  variable,  eqn.  (2) 
subgrid-scale  variable 
dimensional  variable  in  the  Appendix 
three-dimensional ly  Fourier  transformed  variable 
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